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ABSTRACT 


A  method  for  estimating  the  amplitude  and  time  delay 
parameters  for  signals  which  can  be  represented  as  a  sum  of 
a  number  of  scaled  and  delayed  replicas  of  a  known  signal  in 
the  presence  of  non-white  Gaussian  noise  is  presented.  The 
method  is  based  upon  the  principles  of  maximum  likelihood 
estimation  but  at  certain  key  points  simplifying  assumptions 
are  made  that  result  in  a  computationally  practical 
estimation  scheme.  The  results  are  compared  with  other 
related  results  in  the  literature. 
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I.  INTRODUCTION 


A  problem  of  practical  concern  that  has-  received  a  great 
deal  of  attention  is  that  of  time  delay  estimation.  Most  of 
this  past  research  [1],  [2],  [3],  [4],  has  focused  on 
problems  of  the  form 

r:(t)  =  s(t)  +  n1(t) 

r 2 ( t )  -  s(t-D)  +  n2(t) 

Another  more  complicated  problem  of  research  interest  is  one 
in  which  the  model  is 

r  1  ( t )  =  s  ( t ) 

M 

r2(t)  «  n(t)  +  E  A^  s(t-D^) 
i-1 


where  s(t)  and  M  are  known  and  we  desire  to  estimate  the 
scaling  amplitudes  A^,  A2 . . • ,  AM  and  delays  ,  D2,...,DM. 
For  example,  this  could  be  a  problem  where  a  known  probe 
signal' s(t)  excites  a  linear  time  invariant  filter  whose 
output  is  observed  in  the  presence  of  non-white  Gaussian 
noise,  and  where  the  impulse  response  of  the  filter  is  a 
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finite  sum  of  time-shifted  delta  functions  with  different 
amplitudes . 

Closely  related  problems  have  been  studied  and  reported 
in  [5],  [6]  ,[7]  and  [8].  Here,  we  take  a  different  approach 
in  the  derivation  that  results  in  a  practical,  easily 
interpretable,  and  computationally  fast  solution  to  the 
problem.  We  view  this  as  being  a  useful  contribution  to  an 
important  engineering  problem. 

II .  APPROACH/SOLUTION 


First,  we  are  given  that  s(t)  and  n(t)  are  time-limited 
(0<t£T)  sample  functions  of  non-white,  zero  mean  Gaussian 
processes  uncorrelated  with  each  other,  Di:>0,  and  D1  *  0. 
Furthermore,  s(t)  is  completely  known  as  is  the  power 
spectral  density  of  its  process,  G  (f).  We  also  require 
that  the  maximum  difference  between  any  two  delays  is  less 
than  the  total  observation  time.  That  is 

max|D^-Djl  <T  for  i^j 
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Then  we  derive  the  maximum  likelihood  (ML)  estimator  making 
approximations  and  dropping  unnecessary  terms. 

Using  the  results  of  Carter  [2],  we  carr  show  that  the  ML 
estimator  is  to  minimize 


through  the  selection  of  ' s  and  Di's.  The  letter  G 
denotes  power  spectra,  f  is  frequency  and  the  hat  denotes  an 
estimate  based  on  observed  data.  Using  the  fact  that  C ^  is 
the  magnitude  squared  coherence,  letting  L  replace  the 
entire  summation  term  and  switching  the  minimization  to  a 
maximization  by  changing  signs,  we  find  should  maximize 


where  we  have  dropped  the  explicit  notational  reminder  of 
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By  expanding  and  dropping  terms  not  dependent  upon  the  ' s 
and  D^'s  we  obtain  the  new  expression  to  be  maximized  as: 


As  will  be  seen,  this  frequency  domain  formulation  is 
closely  related  to  the  time  domain  formulation  of  Ehrenberg, 
et.  al .  [6],  when  the  appropriate  assumptions  are  made. 


Estimating  the  Amplitudes 


We  observe 
for  estimating 


A 

in  J_  that  G, ,  -  G_  and  test  how  J„  works 
c  1 1  s  s  c 

amplitudes  in  the  following  simple  example: 


r  x  ( t )  =  s  ( t ) 

r 2 ( t )  -  A1  sft-D^  +  n ( t ) 


4 


TM  No.  861057 


Assume  (for  this  example)  that  D^O  and  n(t)  is  white  with 
unit  power.  Then  we  maximize: 

f 4 Afd-P 

“o*  —  ofr 

=  3  Ai  Ria(W-  A? 


A 

where  R^2  is  the  estimated  cross-correlation  function 
between  r1(t)  and  r2(t),  and  Rss  is  the  autocorrelation 
function  of  s( t) . 

Taking  the  partial  derivative  with  respect  to  A^  and  setting 
it  equal  to  zero,  we  get: 

=  3  RqOO'iJAi  ^£s(®)  =  0 

So  tkat 

A|"  R«(°) 


A 

In  this  example,  if  the  noise  were  not  white,  R^2^D1^ 
would  be  the  noise  prewhitened  generalized  cross- 
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correlation  function  evaluated  at  the  delay,  ,  and  which 
we  designate  as  R^2'{D^)/  and  Rgs(0)  would  be  replaced  by 
the  noise  prewhitened  auto-correlation,  R  '(0),  where 

Q  S 


and 


Kss'(°)  - 


<&o 

r 


Gss<f» 


J 


Gnn(f> 


df 


These  results  can  be  generalized  by  taking  the  partial 
derivatives  of  Jr  with  respect  to  A.  and  setting  each  equal 

j 

to  zero. 


^  5 


[s  A)  R„'<  y] 


OO 


3  G^}_ 

Gnn  (0 


df  -  0 

/3 2-j 
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2 

Before  taking  the  partial  derivative  of  |L(f)|  we  identify 

2 

the  critical  portions  of  |L(f)|  depending  on  A j  that  have  a 
non-zero  partial  derivative.  This  yields 


and  lienee 


In  order  to  make  useful  progress  from  here,  we  further 
restrict  our  problem  to  the  case  when  we  can  resolve  the 
delays  relative  to  the  correlation  time  of  s(t).  That  is, 
the  difference  between  any  two  delay  values  must  be  greater 
than  the  correlation  time  of  s(t).  To  help  clarify  this,  we 
determine  the  cross-correlation  between  r1(t)  and  r2(t): 


R12(d)  -  E  [ r 1 ( t )  r2(t+d)] 


M 

-  E  [s(  t)  [  I  Ai  sU-D^d)  +  n(t+d)]] 
i-1 


M 

“  2  Ai  Rss(d-Di) 
i-1 
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So  we  see  that  is  a  sum  of  shifted  replicas  of  Rgs(d) 

with  amplitude  weightings  A..  Now  if  R  (d)=0  for  [d|>d  /2 

1  s  s  o 

we  say  that  values  of  s(t)  are  uncorrelated  for  time 

differences  greater  than  d  ,  and  loosely  refer  to  d  as  the 

c  c 

correlation  time  of  s(t).  Thus  we  can  state  our  restriction 
on  the  differences  in  the  delays  as 


|Di-DJ!>dc  i/j 

With  this  restriction  we  can  make  the  observation  that: 


Rss'  (Dj-V  :  1/3 


Therefore,  we  conclude  that 


^  =  aR,'(D/)-aAjR,>)  =  0 


so  that 


S12'(V 


Aj  - 


Rss'(01 


Thus  the  only  thing  we  need  to  estimate  the  attenuations  is 
the  noise  prewhitened  generalized  cross  correlation  function 
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at  the  correct  time  delays  and  the  noise  prewhitened  signal 
powe  r . 


For  nearly  white  noise 


Rss' 


(0) 


Rss(0>/Gnn 


where  G  is  the  noise  power  spectrum  level.  In  this  case 
nn 


r12'(V  Gn„ 


Rss(0) 
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Estimating  the  Delays 


We  now  obtain  equations  foe  D.'s.  Expand  J_  to  get 

D 


-oo 


-90 


We  note  that  maximizing  JD  is  equivalent  to  maximizing  J^. 
At  this  point  we  derive  the  cross  power  spectrum: 

G12(f)  -  F[R12  (t)] 

M 

=  F [ E ( S ( t )  Z  A.  s(t-Di  +x) } ] 
i-1 
M 

=  Gs$(f)  Z  Ai  exp(-j2nfD. ) 
i-1 

-G$s(f)  L  (  f ) 

Where  F  denotes  taking  the  Fourier  Transform,  and  E  is  the 
expectation . 
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Now  we  focus  attention  on  the  last  two  terms  of  JD>  If 

we  assume  r^(t)  is  a  noise  free  copy  of  s(t)  then  we  do  not 

need  to  estimate  G^(f)  since  it  is  simply  G  (f).  When  we 

get  in  the  neighborhood  of  the  correct  estimate  of  G12  such 
*  *  , 

that  =  GSSL  L  we  see  that  JD  is  approximately  equal 

to 


o O  A 


<re  =i 

-O o  (r, 


An 


We  note  that  this  step  is  an  approximation  since  we  have 
only  time  limited  signal  and  noise  sample  functions.  With 
the  approximation,  the  last  two  terras  cancel.  Hence,  for 
problems  of  practical  interest,  the  function  to  be  maximized 
can  be  written  as 


Even  though  looks  complicated  we  will  now  show  how  this 
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form  of  the  answer  leads  to  a  significant  simplification  in 
the  practical  engineering  solution  of  this  problem. 

The  maximization  of  J„  can  be  rewritten  as 

F 


We  see  that  the  term  to  be  maximized  can  be  interpreted  as  a 
weighted  finite  sum  of  generalized  cross-correlation 
functions  with  noise  prewhitening.  Note  also  that  the 
integration  term  is  a  Fourier  transform  of  a  weighted  cross 
power  spectrum;  once  it  has  been  computed  (e.g.,  via  an  FFT ) 
it  does  not  have  to  be  recomputed  for  other  D^s.  In  terms 
of  the  noise  prewhitened  cross-correlation  function  the 
maximization  is  of 

M  A 

JF  -  I  AiRjjMDj) 

i-1 


by  selection  of  D^'s. 


Now  we  substitute  the  expression  for  derived  earlier  and 


we  have 

Jf 
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Finally,  noting  that  R  '(0)  is  just  a  constant,  we  see  that 

b  b 

maximizing  J?  is  equivalent  to  choosing  the  D^'s  which 
maximize 


The  selected  D^'s  are  also  subject  to  the  constraint  that 

for  i  /  j,  then  D ^  ?  Di  and  | D ^  — Di |  must  be  greater  than  the 

correlation  time  of  s(t).  To  put  it  simply,  we  find  the  M 

2 

highest  peaks  of  [R^2'(d)3  that  are  separated  by  at  least 
the  correlation  time  of  s{t).  If  we  make  use  of  our  earlier 
assumption  that  D1=0 ,  and  introduce  an  additional  constraint 
that  for  i>j  we  can  have  a  meaningful  ordering  of  our 

A.  and  Dj  estimates  without  affecting  the  generality  of  our 
solutions.  Note  that  the  squaring  and  summing  of  the 
generalized  cross-correlation  function  takes  care  of 
negative  peaks  possibly  occuring  at  some  of  the  D^'s. 


IV.  COMPARISON  WITH  OTHER  RESULTS 


Here,  we  make  comparisons  between  our  results  and 
previously  published  works  that  deal  with  problems  closely 
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related,  but  not  identical  to  the  one  we  have  just 
addressed. 

Comparison  to  Ehrenberg,  Ewart,  and  Morris  [6] 

It  is  can  readily  be  shown  that  our  results  and  those  of 
[6]  become  identical  when  we  impose  our  constraint  of 

RSg( Di~Di ) for  i/j 

and  make  their  assumption  of  white  noise  with  two-sided 
power  spectral  density  Nq/2.  So  for  more  restrictive 
assumptions  on  delay  resolvability  but  less  restrictive 
assumptions  on  noise  complexity,  we  obtain  an  easy  to 
implement  processor. 


Comparison  to  Knapp  and  Carter  [3] 


For  a  single  time  delay  modelled  by 


r1(t)  -  s ( t)  +  n^  t) 
r 2 ( t )  -  s(t-D)  +  n2(t) 
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and  under  typical  simplifying  assumptions,  Knapp  and 
Carter[3]  have  shown  that  the  ML  delay  estimator  is  the 
delay  at  which  the  generalized  cross-correlator  with  proper 
prefilters  is  a  maximum.  In  particular,  for  the  single 
delay  model  the  solution  was  shown  to  be 

A 

R12(T)  *  F-1 [ W( f )  G12 ( f ) ] 

where  F-1  denotes  the  inverse  Fourier  transform  and  the 
conjugate  product  of  the  prefilters  of  r^(t)  and  r ^ (t)  is 

W(f)  -  H1(f)  H2* ( f )  -  1  C(f ) 

|G12( f ) (  [l  -  c ( f ) ] 

and  C(f)  is  the  magnitude-squared  coherence  (i.e.,  the 
magnitude  squared  cross-power  spectrum  divided  by  the 
product  of  the  two  autospectra).  Now  for  the  model  here 
with  n^(t)  -  0  and  r^Ct)  -  n(t),  we  find 


C(f ) 


(1 


C(f )  ] 


ss 


(f) 


nn 


(f) 


Also  for  the  model  here 

|Gi2<f)|  =  Gsg(f) 
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So  the  estimator  is  to  find  the  delay  that  causes  the 
maximum  of 

A 

R12(t)  =  F_1  [  1  G12(f)1 
Gnn(fl 

This  is  exactly  the  term  in  our  problem  if  there  were 
only  one  delay  and  if  A^  «  1.  There  are  two  other  points, 
of  interest.  First,  Knapp  and  Carter  [3]  introduced  a 
squaring  of  the  correlation  function  in  their  figure  showing 
a  practical  realization  of  picking  a  maxima;  this,  in  part, 
takes  care  of  negative  attenuation  in  one  receiver  channel. 
Second,  if  we  wanted,  we  could  filter  s(t)  by 
H^f)  =  l/[  Gnn(  f )  ] 1//2  and  ^(t)  the  same  function. 

Whalen  [101  has  derived  the  ML  solution  for  a  simple 
case  of  a  single  attenuation.  In  particular,  given 

r(t)  «  A  s(t)  +  n(t) 


where  n(t)  is  white  noise  the  ML  estimate  of  A  is 
A  ^  gr(t>  s(0  dt_ 

A'  £  S’fOdt  "  fa(o) 


so  we  correlate  r ( t )  with  s(t);  no  noise  prefilters  are 
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required,  since  n(t)  is  white;  then  we  read  off  the 
correlation  function  at  zero  delay  and  normalize  by  the 
signal  power.  If  we  expand  the  solution  above  we  see 


A  _  A  Si  +-  dt 
A  ' 


and  that  by  taking  the  expectation  the  answer  is  unbiased. 

Now  if  n(t)  is  not  white  then  from  r(t)  *  As(t)  +  n(t) 
we  form 


r'(t)  =  As ' ( t )  +  n'  ( t ) 


where  r',  s'  and  n'  are  convolution  outputs  with  a 
prewhitening  filter 


H (  f)  | 


Gnn<f> 
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Then 


a  [J f'W  s'fe)dt  ^  G„,  Ijr  '(i) 

A=  IjPfdt  Rsrf°) 


so  both  s(t)  and  r(t)  are  subjected  to  a  prefilter 


I  H  (  f )  |  2 


1 


'nn 


(f) 


and  thus  the  normalization  takes  into  account  both  the  noise 
and  the  signal  power. 


V.  RE-CAP  OF  STEPS 


1.  Estimate  the  cross  power  spectrum  between  stored  replica 
(reference)  r^(t)  =  s(t)  and  received  signal  ^(t). 

2.  Estimate  the  noise  power  spectrum  Gnn(f)  by  noting 
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g22(£)  '  Gss(£)  lL(£)'2  +  Gnn(£) 


and 


|G12(f)|2 


Gss(f) 


Ggs(f)  |L(f)  |2 


and  computing 


A  A  A 

Gnn(£>  ’  G22<£>  '  I  I G12( £) | 2/Gss( f ) ] 


3.  Compute  the  total  signal  power  R  (0)  and  the  total 

s  s 

noise  prewhitened  signal  power  R__'(0). 

5  b 

4.  Compute  (e.g.,  via  an  FFT )  the  inverse  transform  of  the 
cross  power  spectrum  divided  by  the  noise  power  spectrum  to 
get  the  generalized  cross-correlation  function.  Also 
compute  the  square  of  this  function. 
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5.  Find  the  M  largest  peaks  of  the  squared  generalized 
cross-correlation  function,  which  are  separated  by  at  least 
the  correlation  time  of  s(t).  The  corresponding  delay  values 
at  these  peaks  are  the  M  delay  estimates. 

6.  Determine  the  value  of  the  generalized  cross-correlation 
function  at  each  of  the  M  delay  estimates. 

7.  Adjust  the  generalized  cross-correlation  values  found  in 
step  6  by  dividing  by  the  noise  prewhitened  signal  power. 
These  adjusted  M  correlation  peak  heights  then  become  ML 
estimates  of  the  corresponding  amplitude  scalings.  If  the 
noise  spectrum  is  relatively  fla.t  (i.e.,  nearly  white)  then 
the  estimates  can  be  approximated  by  multiplying  the  values 
found  in  step  6  by  the  noise  power  spectrum  level,  then 
dividing  them  by  the  total  signal  power.  There  should  be 
good  agreement  in  the  two  estimates  for  this  case. 

VI .  SUMMARY 


We  have  derived  a  pratical  method  for  estimating  the 
amplitudes  and  delays  for  a  composite  signal  in  additive, 
non-white  Gaussian  noise,  the  signal  being  composed  of  a  sum 
of  amplitude-scaled  and  time-delayed  replicas  of  a 
completely  known  sample  function  of  a  Gaussian  process  that 
is  uncorrelated  with  the  corrupting  noise.  We  have 


20 


TM  No.  861057 


constrained  the  problem  to  the  case  where  the  composite 
time-delays  are  resolvable  with  respect  to  the  correlation 
time  of  the  basis  signal. 

The  form  of  our  solution  is  easily  interpretable  and 
convenient  to  implement  by  means  of  generalized  correlation 
functions.  We  have  found  that  our  results  are  consistent 
with  related  results  of  previous  work  when  appropriate 
adjustments  of  constraints  and  assumptions  are  made,  while 
our  approach  has  been  based  on  the  principles  of  maximum 
likelihood  estimation,  the  simplifications  made  in  the 
derivation  and  practical  limitations  in  observing  a 
stochastic  signal  make  our  resultant  formulation  approximate 
maximum  likelihood  estimates. 
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